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ABSTRACT 


This  paper  presents  the  development  of  a  finite  difference  form  of  the  unsteady  boundary  layer 
equations  for  twoldimensional  flow.  The  equations  are  written  in  body-conformal  coordinates, 
non-dimensionalized,  transformed  to  a  generalized  curvilinear  coordinate  system,  and  expanded 
to  yield  equations  suitable  for  use  in  a  tridiagonal  solver.  The  boundary  conditions  are  also  presented.  < 


This  work,  done  in  cooperation  with  the  University  of  Toronto,  seeks  the  development  of  a  boun¬ 
dary  layer  code  compatible  with  the  NASA  Ames  ARC2D  Navier-Stokes  code.  The  two  codes  will 
be  used  in  a  study  of  the  Fortified  Navier-Stokes  concept. 


RESUME 


Cette  communication  porte  sur  le  developpement  d’une  forme  aux  differences  finies  des  equa¬ 
tions  des  couches  limites  instables  dans  les  ecoulements  bidimensionnels.  Les  equations  sont  ex- 
primees  en  cooidonnees  conformes  au  corps  etudie,  sans  dimension,  sont  transformees  en  un  systeme 
de  cooraonnees  curvilineaires  generalises  et  sont  developpees  en  equations  applicables  a  un  resolveur 
tridiagonal.  Les  conditions  aux  limites  sont  aussi  presentees. 

Ce  travail,  mene  en  collaboration  avec  l’Universite  de  Toronto,  vise  a  elaborer  un  programme  de 
couche  limite  compatible  avec  le  programme  Navier-Stokes  ARC2D  mis  au  point  par  1’Ames  pour 
la  NASA.  Les  deux  programmes  serviront  k  etudier  le  concept  fortifie  de  Navier-Stokes. 
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NOMENCLATURE 


A,B,C,D 

Coefficients  of  tridiagonal  system  of  equations 

a 

Speed  of  sound 

cp 

Coefficient  of  specific  heat  at  constant  pressure 

H 

Total  enthalpy,  H  =  cp  +  u2H 

l 

Reference  length 

M 

Mach  number,  M-ula 

P 

Pressure 

Pr 

Prandd  number,  Pr  =  pcp/K 

R 

Gas  constant 

Re 

Ratio  of  freestream  Reynolds  and  Mach  numbers.  Re  =  p«. 

Si 

Constant  used  in  Sutherland’s  approximation  for  viscosity 

t 

Time-like  variable  used  to  relax  the  equations 

T 

Absolute  temperature 

u,v 

Body-conformal  velocity  components 

u,v,w 

Contra viariant  velocity  components 

x,y 

Body-conformal  coordinates 

X,Y 

Cartesian  coordinates 

Greek  Symbols 


a  Convective  coefficient,  a  =  p U 

y  Ratio  of  specific  heat  at  constant  pressure  to  specific 

heat  at  constant  volume 

k  Coefficient  of  thermal  conductivity 

p  Viscosity 

p  Density 

x  Shear  stress,  x  =  p(3a  /  By) 

4,T1  Curvilinear  coordinates 


iv 


Oo/Pt» 
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Subscripts 


Oo 

Reference  conditions  in  Sutherland’s  approximation 

Ooo 

Freestream  values 

()/,()* 

Dummy  spatial  indices 

a 

Value  at  a  solid  boundary  (wall) 

(  )wc 

Value  at  the  wake  centerline 

Superscripts 


(  )  Non-dimensionalized  value 

( )"  Dummy  time  index 
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INTRODUCTION 


Background 

Throughout  the  Computational  Fluid  Dynamics  community,  two  approaches  to  the  numerical 
simulation  of  viscous  flow  have  gained  popularity:  the  Interactive  Boundary  Layer  (IBL) 
approach  and  the  Navier-Stokes  approach. 

IBL  methods  solve  separate  inviscid  and  viscid  solutions,  which  influence  each  other  throught  an 
interaction  algorithm.  The  interaction  is  allowed  to  continue  until  continuity  of  the  velocity  and 
temperature  profiles  is  achieved  at  an  interface  boundary,  usually  the  boundary  layer  edge. 
Numerical  implementation  of  the  method  involves  coding  of  separate  solution  algorithms  for  the 
inviscid  and  viscid  regions,  each  based  on  specialized  equations  which  approximate  the  Navier- 
Stokes  equations  under  the  conditions  prevailing  in  that  region.  The  viscous  method  also 
requires  special  adaptations  for  separated  flow  regions,  where  the  boundary  layer  equations  must 
be  solved  in  inverse  mode.  This  results  in  substantial  human  intervention  to  adapt  a  given  code 
to  a  new  class  of  problems. 

Navier-Stokes  methods  solve  the  Navier-Stokes  equations  throughout  the  field.  Their  numerical 
implementation  is  simpler  since  a  single  set  of  equations  applies  to  all  regions  of  the  flow  field, 
but  grid  points  must  be  concentrated  in  the  viscous  regions  to  resolve  the  boundary  layer  flow. 
Since  the  Navier-Stokes  equations  describe  the  physics  for  most  aerodynamic  problems,  minimal 
human  intervention  is  normally  required  to  adapt  a  given  code  to  a  new  class  of  problems. 

It  is  recognized  that  IBL  schemes  require  much  less  computer  work  per  grid  point  and  less 
computer  storage  than  Navier-Stokes  solvers  to  obtain  a  solution.  However,  they  are  more 
specialized  and  are  restricted  to  computation  of  attached  and  mildly  separated  flows.  When 
substantial  flow  separation  exists,  the  less  efficient,  and  more  expensive  Navier-Stokes  methods 
must  be  used. 

Hence,  Navier-Stokes  codes  are  quite  general  in  their  range  of  application  but  they  are  expensive 
to  run.  As  a  result,  most  Navier-Stokes  codes  require  improvements  to  reduce  the  associated 
computational  workload.  On  the  other  hand,  IBL  methods  are  efficient,  but  the  equations  used 
may  not  describe  the  physics  of  the  flow  adequately.  These  codes  require  improvements  to 
produce  more  realistic  simulations. 


Fortified  Navier-Stokes  Approach 

Over  the  past  decade,  NASA  Ames  researchers  used  a  boundary  layer  solution  as  a  forcing 
function  to  the  Thin  Layer  Navier-Stokes  (TLNS)  equations  to  improve  the  convergence 
characteristics  of  their  Navier-Stokes  code.  They  called  this  approach  "Fortified  Navier-Stokes" 
(FNS). 
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Conceptually,  a  FNS  code  is  a  standard  Navier-Stokes  code  which  solves  the  TLNS  equations 
throughout  the  computational  domain.  However,  the  TLNS  equations  have  been  modified  to 
incorporate  approximate  solutions  (e.g.  for  a  boundary  layer  equations)  in  zones  where  they 
deary  apply.  As  a  result,  the  FNS  code  can  be  used  in  a  way  similar  to  an  IBL  code.  An  inviscid 
solution  is  obtained  by  solving  the  TLNS  equations  on  a  coarse  grid.  Coarseness  of  the  grid 
prevents  accurate  resolution  of  the  viscous  terms  within  the  boundary  layer  region.  To  recover 
this  information,  the  boundary  layer  equations  are  solved  using  a  fine  mesh  overlay  near  the 
airfoil  surface,  and  the  solution  is  impressed  on  the  inviscid  solution  through  a  forcing  function 
term  added  to  the  standard  TLNS  equations.  The  advantage  of  a  FNS  code  is  that,  should  the 
boundary  layer  assumption  break  down,  the  forcing  term  can  be  set  to  zero  and  computations  can 
continue  using  the  conventional  Navier-Stokes  solution  on  a  refined  grid.  Steger  and  Van 
Dalsem  developed  several  FNS  code  versions  in  1985^,  1986^  and  1988^,  which  differ  by  the 
form  of  the  boundary  layer  equations  used. 


Previous  Generalized  Boundary  Layer  Work 

Several  boundary  layer  methods  (e.g.  Cebeci^)  employ  coordinate  transformation,  or 
"stretching"  of  the  governing  equations  prior  to  formuladon  of  the  finite  difference  equations. 
The  purpose  is  to  scale  the  growth  of  the  viscous  layer  so  that  the  new  equations,  expressed  in 
terms  of  similarity  variables,  may  be  solved  on  an  approximately  rectangular  grid.  Van  Dalsem 
and  Steger151  took  a  different  approach.  They  developed  a  scheme  to  solve  the  two-dimensional, 
steady  state  boundary  layer  equations  in  body-conformal  coordinates  ( x,y ),  which  amounts  to 
solving  the  equations  in  physical  -  instead  of  similarity  -  variables.  An  adaptive  grid  must  then 
be  used  to  concentrate  a  sufficient  number  of  computational  points  within  the  boundary  layer. 
Solution  of  the  equations  on  a  grid  with  uneven  spacing  results  in  complicated  finite  difference 
equations,  but  this  is  avoided  by  transforming  the  equations  to  a  computational  domain  with 
uniform  grid  spacing  so  that  standard,  unweighted  differences  can  be  used. 

The  scheme  of  Van  Dalsem  and  Steger  is  noteworthy  for  two  reasons.  First,  it  assumes  the 
equations  are  only  weakly  coupled  when  the  x*,,  is  given.  This  allows  them  to  solve  each 
equation  independently  in  a  sequential  scheme.  They  also  chose  their  finite  difference  operators 
to  allow  use  of  efficient  bidiagonal  and  tridiagonal  solvers  in  the  solution  process. 

Second,  their  scheme  avoids  use  of  complex  space  marching  procedures  within  separated  flow 
regions.  Instead,  the  scheme  marches  in  the  general  downstream  direction  while  flow-dependent 
finite  difference  operators  adapt  locally  to  respect  the  parabolic  nature  of  the  boundary  layer 
equations. 

The  earlier  work  of  Steger  and  Van  Dalsem*1^5'  used  a  predictor-corrector  finite  difference 
approximation  to  the  two-dimensional,  steady  state  equations.  They  also  used  three  types  of 
operators  depending  on  the  local  value  of  u/ue.  In  1986'2^  and  1987^',  they  expanded  the 
equations  to  the  three-dimensional  unsteady  form.  They  also  simplified  the  formulation  to  a 
single-step  process  and  only  two  operators  dependent  on  the  sign  of  the  contravariant  velocity  U 
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or  W.  The  time  variable  is  not  used  in  a  true  time  fashion,  since  the  method  still  assumes  that  the 
pressure  distribution  is  fixed  and  given.  Instead,  the  "time-like"  variable  is  used  to  relax  the 
equations.  The  magnitude  of  the  "time  step"  affects  the  convergence  rate. 

In  1983,  Steger  and  Van  Dalsem^  developed  a  completely  different  approach  to  solving  the 
boundary  layer  equations.  The  boundary  layer  equations  were  recast  in  a  form  using  the 
Cartesian  velocity  components  instead  of  the  usual  boundary  layer  components.  They  can  then 
be  solved  on  a  Cartesian  grid,  thus,  increasing  commonality  with  Navier-Stokes  codes. 

In  the  present  work,  the  1986  version  of  the  generalized  boundary  layer  equations  is  used.  It  was 
selected  because  it  is  simpler  to  implement  than  the  !988  formulation  and  will  provide  a 
framework  from  which  to  study  numerical  convergence  issues.  However,  the  equations  are 
scaled  differently  than  in  references  3  and  6  to  increase  their  compatibility  with  Pulliam’s^ 
Navier-Stokes  code,  ARC2D,  with  which  the  current  generalized  boundary  layer  code  will  be 
integrated.  The  equations  are  also  expanded  in  sufficient  detail  for  use  directly  in  a  computer 
program.  Handling  of  the  boundary  conditions  is  also  presented. 
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GOVERNING  EQUATIONS 


Equations  in  Body-Conformal  Coordinates 

We  now  set  to  write  the  two-dimensional  form  of  the  unsteady,  compressible  boundary  layer 
equations  in  generalized  curvilinear  coordinates.  We  start  from  the  dimensional  Cartesian  form 
of  the  equations  as  developed  in  standard  fluid  dynamics  texts^ 


'  du 

—  t  1 1 

du  du 

=  Jte-H 
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du 

t3' 

dx  +V  dy 
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h  ay 

^  ay 
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P  cp 


dT  dT  dT 

-f-  —  +  v  - 

__a_ 
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Ktt- 

+  p 

'  du 

dt  dx  dy 

w  J 

ay 

ayj 

dy 

+  u 


dp_ 

dx 


(Lb) 


ap  3(pu)  3(pv)  _n 

dt  dx  dy 


In  addition,  we  require  two  constitutive  equations:  the  perfect  gas  law,  and  a  relation  linking 
viscos’ty  v  ith  temperature 


p  =  pRT 

(l.d) 

•p 

II 

■p 

(l.e) 

For  air,  the  relation  between  viscosity  and  temperature  can  be  expressed  by  an  interpolation 
formula  based  on  D.M.  Sutherland’s  theory  of  viscosity^8,  p'328^  given  as 


JL_ 

T 

312 

1  o  +  S  i 

Po 

To 

r  +  Sj 
-  - 

where  po  denotes  the  viscosity  at  the  reference  temperature  Tq.  Si  is  a  constant  which,  for  air, 
is  S  i  =  1 10  AT.  This  relation  holds  true  for  subsonic  and  supersonic  flow,  but  not  for  hypersonic 
flow. 
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The  independent  variables  x  and  y  are  the  Cartesian  coordinates;  thus  equation  1  applies  to  flow 
over  a  flat  plate.  However,  by  mapping  the  Cartesian  grid  to  a  body-conformal  grid  (see  fig.l), 
the  functional  form  of  equation  1  remains  valid  as  long  as  the  surface  has  mild  curvature  (e.g. 
the  rounded  leading  ed^  of  an  airfoil,  but  not  around  the  trailing  edge).  In  a  body-conformal 
grid,  the  x  coordinate  is  defined  as  the  arc  length  along  a  solid  surface  (or  along  the  dividing 
streamline  of  an  airfoil  wake),  while  the  y  coordinate  is  the  distance  normal  to  the  body  surface 
(or  wake-dividing  streamline)  and  is  defined  positive  in  the  outward  direction.  The  Cartesian  to 
body-conformal  mapping  of  the  grid  assumes  the  range  of  the  y  coordinates  is  much  smaller  than 
the  range  of  the  x  coordinates.  The  difference  in  arc  length  between  the  lines  defining  the  body 
surface  (y  =  0)  and  the  line  defining  the  outer  edge  of  the  Cartesian  grid  (y  =  ymax)  is  negligible. 


Cartesian  Body-Conformal 


Figure  1.  Cartesian  and  Body-Conformal  Coordinates 

In  their  formulation,  Steger  and  Van  Dalsem'2"6'  use  the  total  enthalpy  form  of  the  energy 
equation  (l.b).  This  form,  derived  in  Appendix  A,  is 


dH 

* 

dH  dH 

_d_ 

JL 

f  dH  Pr-  1 

d(u)2 

dt  +U 

dx  V  dy 

* 

dy 

Pr 

dy  2 

dy 

_ 

(Lb) 


Non-Dimensionalization  of  the  Equations 

To  increase  compatibility  between  the  generalized  boundary  layer  and  Navier-Stokes  codes, 
equation  1  is  non-dimensicualized  with  the  scaling  used  by  Pulliam'7'  in  ARC2D.  Density, 
viscosity  and  temperature  are  scaled  by  their  freestream  values  as 


T  = 


T 


T 

•  CtO 
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Velocity,  enthalpy  and  pressure  are  scaled  by  combinations  of  the  freestream  speed  of  sound 
and  density  as  follows 


;  p  = 


p 

2 

Poo^oo 


Finally,  the  spacial  coordinates  are  scaled  by  a  reference  length  (usually  the  airfoil  chord),  /, 
and  time  is  scaled  by  the  same  quantity  divided  by  the  freestream  speed  of  sound,  i.e. 


x  = 


£ 

/ 


t 

l!a„ 


Using  the  above  scaling,  equations  (l.a-c)  retain  the  same  functional  form  except  that  the 
viscous  terms  of  equations  (l.a,b)  ate  now  divided  by  the  ratio  of  Reynolds  number  to  Mach 
number,  defined  as 


Re  = 


Pflo^  00/ 

Hoc 


Note  that  the  formal  definition  of  the  Reynolds  number  uses  the  freestream  velocity  instead 
of  the  speed  of  sound  a„.  The  standard  Reynolds  number  is  obtained  by  multiplying  the  above 
non-dimensional  Re  by  the  Mach  number  M In  the  remainder  of  this  document,  we  drop  the  ~ 
for  simplicity.  The  non-dimensional  form  of  equations  (l.a-c)  is  then 


■N 

Bu  du  Bu 

05 

1 

1 

Bu 

dt  U  Bx  V  By 

Bx  Re  By 

^ay 

BH  ,  BH 

1  d 

BH  Pr- 1 

a(«)2] 

Bt 

Bx  V  By 

Re  By 

Pr 

By  2 

By  J 

ap  8(p«)  9(py)  n 

dt  dx  By 

The  equation  of  state  for  a  perfect  gas  ( 1  .d)  becomes 


(2.a) 

(2.b) 

(2.c) 
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(2.d) 


where  T  =  (y-  1) 


H  - 


(u2  +  v2) 


and  y  is  the  ratio  of  specific  heat  ( y  =  1 .4  for  air). 


Because  of  its  empirical  nature,  equation  (l.e)  cannot  be  expressed  in  terms  of  scaled  variables 
only.  However,  if  the  freestream  temperature  is  given,  the  following  form  is  obtained 


_ nr  3/2 


\i  =  T 


l  +  Si/T,, 


r  +  Sj/r. 


(2.e) 


Generalized  Curvilinear  Coordinates 

A  general  curvilinear  transformation  is  used  to  map  the  body-conformal  coordinates  ( x,y )  to 
computational  coordinates  (^,rj),  such  that  spacing  in  the  computational  domain  is  uniform  and 
of  unit  length  (see  fig.2).  Equations  (2.a-c)  are  then  transformed  to  general  curvilinear 
coordinates 


$  =  $(*) 

T1  =  T](x,y) 

and  solved  on  the  computational  domain  using  standard  unweighted  finite  differences.  The 
advantage  of  this  transformation  is  that  a  single  computational  code  can  be  used  to  solve  the  flow 
about  a  wide  variety  of  physical  geometries. 

There  normally  exists  a  one-to-one  correspondence  between  the  body  conformal  and 
computational  spaces.  A  notable  exception  is  the  wake  cut  when  a  C-grid  geometry  is  applied  to 
airfoils  (see  fig.3).  Since  the  body-conformal  space  is  obtained  by  "cutting"  physical  space 
along  the  wake  and  "unwrapping"  it  to  align  the  wake  centerline  and  body  surface  along  the 
body-conformal  x  axis,  physical  points  along  the  wake,  including  the  airfoil  trailing  edge, 
correspond  to  two  boundary  points  in  the  computational  domain.  It  is  important  to  choose  the 
wake  points  carefully.  They  should  correspond  closely  to  the  physical  wake  centerline. 
Application  of  the  boundary  conditions  for  the  current  boundary  layer  algorithm  is  then  greatly 
reduced,  as  will  be  discussed  later. 
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Body-Conformal 


Computational 


Figure  2.  Body- Conformal  and  Computational  Grids 


Cartesian  C-Grid 


Computational  Grid 


Figure  3.  Mapping  of  C-Grid  to  Computational  Domain 


The  chain  rule  is  used  to  express  the  body-conformal  derivatives  ( dx  and  dy )  of  equations  (2.a-c) 
in  terms  of  the  curvilinear  derivatives 


dx  —  +  'Hx^ri 

dy  — 

Transformation  of  equations  (2.a-c)  using  these  derivatives  yields 
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ut  +  Uu §  +  =  -4* p 5  + 


pfff.  +  w/t  +  wj  =-^ny 


JL 

Pr 


Pr  -  1 


Tly^ti  +  - j - V“  *1 


Pr  +  ^(PM)4  +  tL(Pw)t,  +  r\y(pv\  =  0 


(3.a) 


(3.b) 


(3.c) 


where  U  and  V  arc  the  contravariant  velocities  defined  as 


U  =  ^u 


V  =  ^yu  +  Tly  v 


The  equation  of  state  (2.d)  and  the  Sutherland  formula  (2.e)  are  not  affected  by  the 
transformation. 

Equations  (3.a-c)  are  similar  to  those  of  references  3  and  6  with  the  exception  that  scaling  of  the 
pressure  and  viscous  terms  is  changed.  In  the  next  section,  we  apply  the  1986  version  of  the 
Steger  and  Van  Dalsem  finite  difference  scheme  to  the  current  equations  and  obtain  the  finite 
difference  equations  required  for  coding. 
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NUMERICAL  FORMULATION 


Steger  and  Van  Dalsem^^  developed  a  time-like  algorithm  to  solve  the  boundary  layer 
equations.  Their  method  assumes  that  pressure  is  known.  This  allows  decoupling  of  the 
equations  which  are  then  solved  sequentially  at  each  time  step  or  cycle.  Furthermore,  the 
algorithm  was  selected  so  that  scalar  tridiagonal  systems  of  equations  are  solved  for  each 
equation.  Time  derivatives  are  first-order  accurate  while  spacial  derivatives  in  the  momentum 
and  energy  equations  are  approximated  by  centered  differences  in  the  q  direction  and  by  flow- 
dependent  upwind  differences  in  the  4  direction.  At  any  £,  station,  the  upwind  direction  is 
determined  from  the  sign  of  the  coefficient  p U.  A  backward  difference  is  used  if  p U  is  positive 
and  a  forward  difference  is  used  if  p U  is  negative. 

The  term  (pv)^  of  the  continuity  equation  is  integrated  in  the  q  direction  using  the 
trapezoidal  rule.  The  £  and  q  derivatives  for  the  (p u)  terms  are  approximated  by  second-order 
accurate  centered  differences. 

Details  of  the  algorithm  are  now  presented  for  each  equation,  separately,  using  conventional 
operators  defined  in  terms  of  the  shift  operator  E,  e.g.E^Uj  =  uJ±\ 


and  the  second-order  accurate  upwind  operator 


a-  lal 
2 


where  a  is  the  convective  coefficient  p U. 


Finally,  we  use  the  convention  that  space  or  time  indices  are  written  only  if  changed,  i.e. 

,,  _  ,.n  ,,rt+l  _  ,.n+l  ~tr 

u  =  Ujj t,  u  =  u j'k  » etc. 
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Momentum  Equation 

General  Form 

Within  the  main  field  (i.e.  away  from  boundaries),  the  finite  difference  form  of  equation  (3.b)  is 
given  by 


pIv-1  +  =-5,8y>  +  ^A(nA»ri 


To  be  solved  numerically,  the  above  equation  must  be  expanded  in  terms  of  the  numerical 
values,  which,  after  rearranging,  leads  to  the  following  form 


pV  i}y 
2  Re 


(4T1>)  +  0^1  A-t 


n+l 

Wjk-1 


+ 


_e_ 

At 


+  p 


U+  l(/l 

3 

V-  It/I 

2, 

h2k_ 

'  (nVjt-i +  +  (^A+i  ] 

2 

-  * 

2 

2 

k.  a 

2  H 

Re 

2 

„  y  _ 

,  n+l 


pV 

2  Re 


(M%)  + 


,.n  +  I 

M*+ 1 


(4.a) 


=  ~kx 


-Pi- \ 

U+  It/I 

-4wy-i  +  M;-2 

r 

§ 

1 

V _ 

4«;+l  -  Uj+  2 

2 

~  P 

2 

a  a 

2 

a  a 

P 

2  J 

2 

a  - 

Equation  (4.a)  is  solved  implicitly  at  each  £,  station.  Thus,  a  tridiagonal  set  of  equations  is 
obtained  by  solving  this  equation  for  the  points  Tj  =  2  to  Tj  =  r\max-\  ■ 

Use  of  equation  (4.a)  near  the  ^  boundaries  is  unacceptable  numerically.  For  example,  at  the 
£=  1  boundary,  the  variables  uj-\ ,  uj-2  and  pj-\  reference  values  outside  of  the  numerical 
arrays  U()  and  P().  This  introduces  errors  and  usually  leads  to  numerical  instability.  At  this 
point,  it  is  useful  to  consider  the  treatment  of  equation  (4.a)  near  the  boundaries  of  the 
computational  domain. 


S,  Boundary  Condition 

The  choice  of  boundary  conditions  depends  on  the  flow  characteristics  of  the  physical  domain 
and  the  choice  of  grid  topology.  The  C-grid  used  in  the  present  work  extends  far  enough  into  the 
wake  to  ensure  that  outflow  conditions  exist  at  the  %  boundaries.  This  means  that  forward 
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differencing  of  the  pU  term  is  used  at  £  =  1  and  backward  differencing  is  used 
Nevertheless,  the  numerics  require  that  the  second-order  accurate  upwind  operator 


y^y^max- 1 

and 

£=  1: 

a8t  - 

^  =  2: 

a8^  = 

a5^  = 

^max: 

a8^  = 

at  ^  =  4max  • 

be  modified 


The  above  modifications  to  the  streamwise  convective  derivatives  ensure  that  only  points  within 
the  computational  domain  are  used.  Hence,  only  first-order  differences  are  used  to  approximate 
inflow  derivatives  at  the  %  =  2  and  E,  =  stations.  Dropping  the  backward  difference  at  the 

£,  =  1  station  and  the  forward  difference  at  the  £,  =  £max  may  not  be  physically  acceptable,  but  it 
is  necessary  numerically.  The  onus  is  on  the  user  to  make  sure  that  outflow  conditions  prevail  at 
these  stations. 

To  test  coding  of  the  algorithm,  other  grid  types  may  be  used  and  inflow  conditions  may  exist  at 
one  of  the  ^  boundaries.  Such  is  the  case  for  the  computation  of  the  flow  over  a  flat  plate.  In 
this  case,  upstream  data  must  be  supplied  (e.g.  for  stations  £  =  1,2  if  pU  is  positive). 

Treatment  of  the  pressure  terms  at  the  2;  boundaries  simply  consists  of  replacing  the  centered 
difference  by  a  forward  difference  at  ^  =  1  and  by  a  backward  difference  at  £  =  £max.  Since  the 
pressure  terms  are  known  a  priori,  they  do  not  affect  the  tridiagonal  nature  of  the  algorithm  and 
second-order  accurate  differences  are  used  as  follows 


j 
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ti  Boundary  Condition 

Two  distinct  boundary  conditions  are  applied  along  the  T|  =  1  line.  The  no-slip  condition, 
(u  i  =  0),  is  used  for  points  on  the  surface  of  the  airfoil.  Thus,  at  the  q  =  2  point,  equation  (4.a)  is 
expressed  as 


A2mi  +B  2m2  +  C2M3  =  D2 

where  A2,B2,  C2  are  the  coefficients  of  ut±\,un+1  andw*+j  respectively  in  equation  (4.a).  D2 
represents  the  right  hand  side  of  the  equation.  With  u  1  set  equal  to  zero,  the  first  term  simply 
disappears  from  the  equation. 

In  the  wake,  the  boundary  condition  uy  I  {  =  0  is  used.  This  derivative  may  be  expressed  as  a 
second-order  accurate  forward  difference  in  q,  i.e. 


Vi=An 


3  ~E 


+1 


un+l  =0 


A  value  for  u  1  in  terms  of  u2  and  M3  is  obtained  from  the  above  equation.  It  is  substituted  in  the 
coefficient  form  of  equation  (4.a)  at  q  =  2,  which  yields 


"  4A2 

r  a*} 

3  +B2 

M2  + 

- - > 

'  1  ^ 

1 

u 

_ / 

U3  -D2  (4a,  1) 


Locally,  the  tridiagonal  form  of  the  system  of  equations  is  preserved. 

At  the  upper  q  boundary  (q  =  qmax),  the  standard  boundary  layer  edge  condition  uy\g  is  used. 

The  same  logic  used  to  derive  the  wake  centerline  boundary  condition  yields  the  following 
coefficient  form  of  equation  (4.a) 


AKm,x-i  ~ 


CKm.x-l 


Ur  -1  + 

^  mix  1 


Br  _i  + 

mix  1 


*CK  _j 

^  mix  1 


uKm„- 1  =  Dr  (4.a,2) 


Inverse  Formulation 

Near  and  in  regions  of  reverse  flow,  the  pressure  gradient  term  of  equation  (4.a)  must  be 
replaced  by  an  inverse  forcing  function.  Steger  and  Van  Dalsem  use  the  wall  shear  stress,  xw,  on 
the  airfoil  surface  and  the  wake  centerline  velocity,  u„c  in  the  wake.  These  expressions  are 
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obtained  by  applying  the  momentum  equation  at  the  airfoil  surface,  which  yields  the  following 
relation 


W 


_2_ 

Re 


Hi  +  H2  “2  "Ml 

2  yi-y\ 


-Xw 


(5.  a) 


Note  that  U\  =0  and  that  xw  is  scaled  by 
relation  is 


/ _ 

(|i~0  ' 


At  the  wake  centerline,  the  equivalent 


kxP^  = 


Hi 

Re 


i]  pUu^ 

J  T1  i 


wc 


(5.b) 


_2_ 

Re 


Hi  +  M-2  M2-«i 

- xw 

2  ya-yi 

y2-),i 


wc 


where  u  j  =  uwc,  t*,  =  0,  and  U I  wc  = 

The  presence  of  a  1*2  term  the  inverse  forcing  functions  results  in  an  augmented  tridiagonal 
system  of  equations.  Steger  and  Van  Dalsem'4'  developed  an  efficient  lower-upper 
decomposition  scheme  to  solve  this  particular  system  of  equations. 


Energy  Equation 

General  Form 

The  finite  difference  scheme  used  for  the  momentum  equation  is  now  applied  to  equation  (3.b), 
the  energy  equation.  Using  the  finite  difference  operators  defined  previously,  we  obtain 


p[v,Hn+{  +  USt>Hn+]  +  V5^Wn+1j  = 


JL 

Pr 


nA  h 


rt  +  1 


Pr- 1  ~ 


TlA(“n+1)2 
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Expanding  this  equation  in  terms  of  the  field  values,  and  rearranging  to  the  tridiagonal  (in  q) 
coefficient  form  ,  we  obtain 


pV  Tiy 
'  2  2  Re 


Pr 


Pr 


k- 1 


jjrt+1 

™k- 1 


A/ 


C/  +  I  (/I 


~P 


U  -  It/I 


2  3_ 

2  2  Re 


Pr 


+  2 


Pr 


Pr 


k+\ 


Hn 


py  By 

2  2  Re 


PBy 


Pr 


Pr 


k+ 1 


f/rt+1 

w*+l 


(4.b) 


By 


2  Re 


Pr  -  1 


PBy 


Pr 


*+l 


PBy 


Pr 


[«Li  -w2]  - 


PBy 


Pr 


k- 1 


1% 


Pr 


2  2 
U  -Uk-\ 


The  solution  scheme  for  equation  (4.b)  is  the  same  as  that  for  equation  (4.a).  An  implicit  system 
of  equations  in  q  is  solved  at  each  ^  station.  Also,  the  restriction  on  the  numerical  values  of 
Hj±  1  and  Hj± 2  near  the  b,  boundaries  are  the  same  as  for  the  momentum  equation. 


Boundary  Conditions 

The  explanation  given  for  the  £,  boundary  conditions  of  the  momentum  equation  also  holds  true 
for  the  energy  equation.  Two  boundary  conditions  may  be  applied  along  the  q  =  1  boundary. 
The  wall  and  wake  centerline  may  be  specified  ( [Tw  and  Twc  given)  or  adiabatic  conditions  may 
be  used  (  Ty\w  and  Ty  I  wc  set  to  zero).  In  the  present  work,  the  adiabatic  condition  is  used 
throughout  the  domain  which  leads  to  the  condition 


Hy  I  w  —  Hy  I  WC 


=  0 


or. 


r 

—a2+b2 

k  - 


H  2  + 


-D2 


(4b,  1) 
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Continuity  Equation 

The  continuity  equation  (3.c)  is  rearranged  to  isolate  the  (pv)  term  on  the  left-hand  side  of  the 
equal  sign  as  follows 


(pvOn  =- 


^x(pU)^+T\x(pu\]+pt 

_ 


The  right-hand  side  is  then  a  function  of  the  variables  p  and  u  alone,  which  are  presumed  known. 
The  trapezoidal  rule  is  used  to  integrate  (pv)  in  the  q  direction  starting  from  the  q  =  1  line  where 
vw  is  known  to  be  zero  on  the  airfoil,  due  to  the  no-slip  condition,  and  vwc  is  assumed  to  be  zero 
due  to  the  symmetry  condition  in  the  wake,  condition  which  should  be  satisfied  if  the  wake  cut 
location  is  selected  carefully.  This  leads  to  the  folowing  finite  difference  equation 

4*8^ (pt*)  4-  q^S^Iptt)  +  Vfp  "+1 


Vn(pvr+I=- 


1  +E 


-l 


Expanding  the  above  equation  in  terms  of  field  variables,  we  obtain  an  explicit  equation  for  (pv) 


(pv)-(pv)ifc_] 


(4.c) 


2  L 


(piO#  -  (pu)fr1 


(p^)!?+i  "  (P“)k-i 

n**-t 


rt+1 


f  r 


(Pu)j+l,k-l  ~  (Pu)]-\,k-\ 


P"+1-P 


A  t 


/r\y 


{pu)n+'-(pu)"kZ2 


n+1 


P*-l  ~Pk-l 


At 


/TW 1 


Centered  differences  in  t,  and  q  are  used  throughout  the  field  to  approximate  the  (pu) 
derivatives.  At  the  boundaries,  centered  derivatives  are  replaced  by  second-order  accurate 
forward  or  backward  differences  to  respect  the  constraints  imposed  by  the  numerical  method,  i.e. 
that  only  points  within  the  computational  domain  be  used.  This  includes  the  following  cases: 


[1] 


At  ^  =  1  '•  Replace  8^(pu)  by  A^ 


3  -  E^' 
2 


(P‘0 
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[2] 

[3] 

[4] 


At  ^  =  4ma* :  Replace  84  (pw)  by 
At  T|  =  2:  Replace  by 

At  rj  =  Tjmax:  Replace  5n(p u)  by  Vn 


3  -El1 


3  -E 


+1  1 


3  -e: 


(p«) 


(pM)*-l 


(PM) 


J 


Application  of  equation  (4.c)  to  points  from  r|  =  2  to  Tj  =  Timax  yields  a  bidiagonal  system  of 
equations  which  is  solved  using  the  standard  tridiagonal  solver.  The  normal  velocity  distribution 
is  finally  obtained  by  dividing  the  result  by  the  density. 


Iterative  Scheme 

The  iterative  scheme  used  in  the  current  implementation  is  identical  to  that  used  by  Steger  and 
Van  Dalsem  in  their  1986-1987  version.  It  is  repeated  here  for  completeness. 

[  1  ]  Solve  the  momentum  equation  (4.a),  throughout  the  computational  domain,  to  yield  the  u 
field  at  time  n+1. 

[2]  Using  the  updated  values  of  u  (which  also  affect  U  and  V),  solve  the  energy  equation 
(4.b)  to  obtain  the  H  field  at  time  n+1. 

13]  Use  the  u”+1  and  Hn+l  values  with  the  equation  of  state  (2.d)  to  obtain  p"+1 . 

[4]  Integrate  the  continuity  equation  (4.c)  to  get  vrt+1  which  is  also  used  to  update  the 
contravariant  velocity  V. 

[5]  Update  the  viscosity  field  p  using  Sutherland’s  formula  (2.e).  For  turbulent  flow,  the 
Baldwin-Lomax  model^9'  is  also  used  to  update  the  values  of  the  apparent  viscosity  p, 
due  to  the  Reynolds  stresses. 

References  3  and  6  are  not  specific  about  the  method  used  to  sweep  the  computational  field.  For 
example,  it  is  not  clear  whether  the  flow  field  is  swept  in  the  general  downstream  direction  from 
the  forward  stagnation  point  or  if  the  sweeps  are  done  from  one  end  of  the  computational  field  to 
the  other  irrespective  of  the  general  flow  direction.  The  equations  developed  in  this  paper  will 
be  used  to  investigate  this. 
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FUTURE  WORK 


The  finite  difference  equations  approximating  the  compressible  boundary  layer  equations  in  two 
dimensions  have  been  developed  in  sufficient  detail  for  use  directly  in  a  computer  code. 
Development  of  a  research  code  to  study  the  convergence  characteristics  of  these  equations  is 
underway.  The  iterative  scheme  outlined  in  section  3.4  will  be  used  with  various  techniques  to 
sweep  the  computational  field. 

The  effect  of  the  time-like  variable  used  to  relax  the  equations,  as  well  as  the  treatment  of  the 
convective  terms  (p U  terms)  in  the  momentum  and  energy  equations,  are  not  well  understood. 
The  new  code  will  be  used  to  study  the  effect  of  these  terms  on  the  overall  convergence 
characteristics  of  the  algorithm. 
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APPENDIX  A 


Enthalpy  Form  of  the  Energy  Equation 

A  rigorous  uevelopment  of  the  enthalpy  form  of  the  energy  equation  is  done  from  the  full 
Navier-Stokes  equations.  However,  in  the  context  of  the  boundary  layer  approximation  to  the 
enthalpy  equation,  the  same  result  is  achieved  by  using  equations  (l.a,b),  which  are  already 
subject  to  the  boundary  layer  assumption. 

Our  development  starts  with  the  addition  of  equation  (l.b)  and  equation  (l.a),  the  latter  being 
multiplied  by  the  velocity  u  as  follows 


bt 

Pcp  a, 


ar  bt 

U  Bx  +v  ay 


a  bt 
a7  Ka7 


+  pw  aT  +  “a7  +  va^  = 


Bu  Bp  Bp  Bp  ,  ?  Bu 

—  +-^-4 ■U~-U^r~  +  U-r-  p.  ~r 

By  Bt  Bx  Bx  By  By 


.  .  ™  ...  Bu  Bu  .  Bu 

We  first  consider  the  left-hand  side  of  the  above  equation.  Terms  like  u  u  —  anil  u  — 


may  be  expressed  as 


B(u2/ 2)  B(u2I2) 


and  u 


B{uz!  2) 


respectively.  Using  the  perfect  gas 


relation  h  =cpT,  the  constant  cp  may  be  brought  into  the  T  derivatives  and  replaced  by  h.  Total 
enthalpy  is  defined  as  H  =  h  +  (u2  +  v2  )/2  ,  but  the  v2  term  is  neglected  under  the  ooundary 
layer  assumption.  Hence  the  left-hand  side  of  the  equation  becomes 


BH  BH  L  BH 

— +  w  — +  v  — 

Bt  Bx  By 


The  right-hand  side  of  the  equation  is  now  considered.  The  u  terms  cancel  out.  The  k 

variable  is  replaced  by  its  Prandtl  number  definition  k  =  \icp/Pr.  We  also  make  use  of  the 
relation 
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a 

du 

__a_ 

d(u2n) 

_  1 1 

du 

M  ay 

^  ay 

<  J 

dy 

dy 

L  J 

M- 

ay 

The  n 


du 

dy 


terms  cancel  out. 


side  of  the  equation  becomes 


Combining  the  terms  within  the  y  derivatives,  the  right-hand 


d_ 

dy 


P-Cp  dT  d(u2/l)}  dp 

Pr  dy  H  dy  dt 

* 


The  constant  cp  is  brought  within  the  T  derivative  and  replaced  by  the  enthalpy.  We  also  add 
and  subtract  the  term  ~~  —  within  the  y  derivative.  Rearranging,  our  equation  (L.H.S.  = 

R.H.S.)  becomes 


dH  dH  dH 

In  |  ii 

_d_ 

dH  ]  Pr  -  1  9(w2) 

dt  dx  dy 

fc  4 

dy 

Pr 

9y  2  dy 

w  4 
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dp 

dt 


Finally,  in  the  current  numerical  scheme,  the  pressure  field  is  used  as  the  forcing  function  to  the 
momentum  equation.  It  is  assumed  known  and  constant,  hence  the  term  is  dropped. 
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